function [bz,az]=bilinear_m(b,a,Fs)
    M=length(b)-1;
    N=length(a)-1;
    L=N-M;
    if L>=0
        b=[zeros(1,L),b];
    else
        a=[zeros(1,-L),a];
    end
    Nz=2*Fs*[1,-1];
    Dz=[1,1];
    N=max(N,M);
    bz=0;
    az=0;
    for k=0:N
        pld=1;
        pln=1;
        for m=0:k-1
            pld=conv(pld,Dz);
        end
        for m=0:N-k-1
            pln=conv(pln,Nz);
        end
        bz=bz+b(k+1)*conv(pln,pld);
        az=az+a(k+1)*conv(pln,pld);
    end
    az=az/az(1);
    bz=bz/bz(1);
end